This is a reanalysis of the main C-peptide result in “Metabolic and immune effects of immunotherapy with proinsulin peptide in human new-onset type 1 diabetes” by Ali et al., Science Translational Medicine 9, 2017-08-09. The main challenges for the primary C-peptide analyses are to properly analyze serial (longitudinal) data within patient, and to obtain a reliable statistical assessment of the evidence for differences among the three treatments.

The R statistical computing system is used for the analysis. See the chapter in BBR on serial data for more details about the statistical analysis used here.

Preliminaries

require(Hmisc)
knitrSet(lang='markdown', w=6, h=4)
require(rms)
options(prType='html')   # for anova and printing fit results
d <- csv.get('dat.txt', header=TRUE, sep='')
d <- upData(d, baseline=baseline/100,
            m3=m3/100, m6=m6/100, m9=m9/100, m12=m12/100,
            tx=c(rep('hf', 9), rep('lf', 10), rep('placebo', 8)))
Input object size:   1968 bytes;     6 variables     27 observations
Modified variable   baseline
Modified variable   m3
Modified variable   m6
Modified variable   m9
Modified variable   m12
Added variable      tx
New object size:    2888 bytes; 7 variables 27 observations
base <- subset(d, select=c(ptno, baseline))

# Reshape from wide to tall and thin dataset
z <- reshape(d, idvar='ptno', varying=c('baseline', 'm3', 'm6', 'm9', 'm12'), v.names='cpeptide',
             timevar='week', times=c(0, 3, 6, 9, 12), direction='long')

i <- with(z, order(ptno, week))
z <- z[i, ]

The first step in any analysis is to plot raw data, especially when the number of subjects is not huge. For serial data we use a “spaghetti plot” of each patient’s time-response profile, starting with a graph on the original scale for C-peptide. Besides never showing such a display, another problem with the paper is its use of normalization to the baseline value, which covers up some possible important features of the data such as whether patients starting very high or very low had different trajectories over time. Always plot the rawest form of the data, with no normalization, and seldom use normalization in the analyses either. Treatments are distinguished by colors in the plot below.

ggplot(z, aes(x=week, y=cpeptide, group=factor(ptno), color=tx)) + geom_line()

In the above graph, the original scale results in compression of the curves making it difficult to see trends, and gives too much influence to a couple of patients with very high C-peptide levels. Let’s redraw the graph using log C-peptide.

ggplot(z, aes(x=week, y=log10(cpeptide), group=factor(ptno), color=tx)) + geom_line()

The influence of the high values is much less, the the distribution of the values on the y-axis direction informally looks suitable for an analysis that assumes normality. But the patient-to-patient (biologic) variability is so great that only a much larger sample size would have made it possible to draw any conclusions from the C-peptide data.

Next merge baseline C-peptide with all the follow-up values, preserving a tall and thin data format.

z <- subset(z, week > 0)
z <- merge(base, z, by='ptno')

Formal Analysis Using A Serial Data Model

The simplest generalization of the t-test or ANOVA to serial data is generalized least squares. We use the AR1 correlation structure, which has been found to fit serial data very often. This specifies an exponentially declining correlation between two measurements on the same patient, as a function of the time between the measurements.

For such a small number of patients, the only hope to analyzing serial data is to assume a smooth time-response mean profile. Here we assume the curves are quadratic (which contains linear trends as a special case). We model treatment effect by interacting the three treatments with the linear and quadratic terms, and we adjust for log baseline (we are analyzing log follow-up C-peptide). The test for treatment is a “chunk” test (a multiple degree of freedom test) of all the terms in the model involving treatment. This tests for a treatment main effect or a treatment by time interaction effect. This could also be interpreted, more naturally, as testing whether there is a difference between any two treatments at any time point.

require(nlme)
dd <- datadist(z); options(datadist='dd')
f <- Gls(log10(cpeptide) ~ pol(week, 2) * tx + log10(baseline),
         correlation=corCAR1(form = ~ week | ptno), data=z)
f
Generalized Least Squares Fit by REML 
 Gls(model = log10(cpeptide) ~ pol(week, 2) * tx + log10(baseline), 
     data = z, correlation = corCAR1(form = ~week | ptno))
Obs 99 Log-restricted-likelihood -24.04
Clusters 27 Model d.f. 9
g 0.322 σ 0.2976
d.f. 89
&beta^ S.E. t Pr(>|t|
Intercept  -0.1246  0.1990 -0.63 0.5329
week  -0.0119  0.0564 -0.21 0.8334
week2   0.0000  0.0036 0.00 0.9963
tx=lf  -0.0647  0.2739 -0.24 0.8138
tx=placebo  -0.2684  0.2813 -0.95 0.3426
baseline   1.0284  0.1836 5.60 <0.0001
week × tx=lf   0.0518  0.0772 0.67 0.5042
week2 × tx=lf  -0.0033  0.0050 -0.66 0.5121
week × tx=placebo   0.0278  0.0786 0.35 0.7246
week2 × tx=placebo  -0.0011  0.0051 -0.22 0.8293 
 Correlation Structure: Continuous AR(1)
  Formula: ~week | ptno 
  Parameter estimate(s):
       Phi 
 0.8757248 
 
anova(f)
Wald Statistics for log10(cpeptide)
χ2 d.f. P
week (Factor+Higher Order Factors) 2.25 6 0.8949
 All Interactions 0.82 4 0.9360
 Nonlinear (Factor+Higher Order Factors) 1.01 3 0.7982
tx (Factor+Higher Order Factors) 4.84 6 0.5646
 All Interactions 0.82 4 0.9360
baseline 31.39 1 <0.0001
week × tx (Factor+Higher Order Factors) 0.82 4 0.9360
 Nonlinear 0.45 2 0.7969
 Nonlinear Interaction : f(A,B) vs. AB 0.45 2 0.7969
TOTAL NONLINEAR 1.01 3 0.7982
TOTAL NONLINEAR + INTERACTION 1.38 5 0.9262
TOTAL 43.59 9 <0.0001

The key test is the overall chunk test for treatment with 6 d.f. The p-value is 0.56, indicating no evidence from these samples for any effect of treatment at any time. Note that the authors also found no overall differences, but then dangerously proceeded to look for an effect at each time, without adjusting for multiple tests. They found one time point where there was evidence for a treatment effect. Note that if one wants to be able to make a statement at any time point, it is important to first do the chunk test for all times combined.

Now lets plot the predicted means (after logging the raw data) as a function of follow-up time, using the fitted quadratic equations, along with 0.95 pointwise confidence bands. These bands can be turned off by clicking in the legend on the right.

wks <- seq(3, 12, length=100)
plotp(Predict(f, week=wks, tx))

Computing Environment[^2]

 R version 3.4.2 (2017-09-28)
 Platform: x86_64-pc-linux-gnu (64-bit)
 Running under: Ubuntu 17.10
 
 Matrix products: default
 BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
 LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
 
 attached base packages:
 [1] stats     graphics  grDevices utils     datasets  methods   base     
 
 other attached packages:
 [1] bindrcpp_0.2    nlme_3.1-131    rms_5.1-2       SparseM_1.77   
 [5] Hmisc_4.0-4     ggplot2_2.2.1   Formula_1.2-2   survival_2.41-3
 [9] lattice_0.20-35
To cite R in publication use:

R Core Team (2017). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.